% AN 169 LINE 3D TOPOLOGY OPITMIZATION CODE BY LIU AND TOVAR (JUL 2013)
function Copy_of_top3d(nelx,nely,nelz,volfrac,penal,rmin)

% 用户自定义迭代参数
maxloop = 200;    % Maximum number of iterations
tolx = 0.01;      % Terminarion criterion
displayflag = 1;  % Display structure flag

% 用户自定义材料属性
E0 = 1;           % 实心材料杨氏模量，最大刚度
Emin = 1e-9;      % 空洞材料杨氏模量，最小刚度
nu = 0.3;         % 泊松比：材料在一个方向受力时，垂直方向的相对形变

% 用户自定义载荷节点自由度
[il,jl,kl] = meshgrid(nelx, 0, 0:nelz);                 % 载荷节点网格坐标：结构右端下边缘垂直载荷（x为右端点，y为0，z遍历所有层）
loadnid = kl*(nelx+1)*(nely+1)+il*(nely+1)+(nely+1-jl); % 节点全局id
loaddof = 3*loadnid(:) - 1;                             % y方向自由度编号：垂直载荷方向

% 用户自定义固定支撑节点自由度
[iif,jf,kf] = meshgrid(0,0:nely,0:nelz);                  % 固定支撑节点坐标：左端靠墙的一侧节点（x=0,y,z遍历所有层）
fixednid = kf*(nelx+1)*(nely+1)+iif*(nely+1)+(nely+1-jf); % 固定支撑节点全局id
fixeddof = [3*fixednid(:); 3*fixednid(:)-1; 3*fixednid(:)-2]; % 固定支撑节点自由度编号：3个方向自由度

% PREPARE FINITE ELEMENT ANALYSIS
nele = nelx*nely*nelz;                      % 单元总数量
ndof = 3*(nelx+1)*(nely+1)*(nelz+1);        % 自由度总数量：节点数量在每个方向上比单元数量多1，一个节点有3个自由度
F = sparse(loaddof,1,-1,ndof,1);            % 节点载荷向量初始化：ndof行1列，载荷节点对应的自由度编号位置大小为1
U = zeros(ndof,1);                          % 节点位移向量初始化：ndof行1列，初始值为0
freedofs = setdiff(1:ndof,fixeddof);        % 无约束节点自由度编号：所有节点-固定节点（差集）
KE = lk_H8(nu);                             % 单元刚度矩阵：
nodegrd = reshape(1:(nely+1)*(nelx+1),nely+1,nelx+1);                   % xy平面中第一个网格平面的所有节点ID：z=0
nodeids = reshape(nodegrd(1:end-1,1:end-1),nely*nelx,1);                % xy平面内部节点编号：列向量，去除nodegrd最后一行，最后一列
nodeidz = 0:(nely+1)*(nelx+1):(nelz-1)*(nely+1)*(nelx+1);               % z方向上不同层第一个节点编号
nodeids = repmat(nodeids,size(nodeidz))+repmat(nodeidz,size(nodeids));  % 三维空间中所有节点编号
edofVec = 3*nodeids(:)+1;                                               % 列向量：每个节点第一个自由度编号
edofMat = repmat(edofVec,1,24)+ ...                                     % 单元连接矩阵：nele*24，每行表示一个单元的24个自由度编号
    repmat([0 1 2 3*nely + [3 4 5 0 1 2] -3 -2 -1 ...
    3*(nely+1)*(nelx+1)+[0 1 2 3*nely + [3 4 5 0 1 2] -3 -2 -1]],nele,1);
iK = reshape(kron(edofMat,ones(24,1))',24*24*nele,1);                   % 行索引：单元刚度矩阵KE的每一行填充到全局刚度矩阵K中的位置
jK = reshape(kron(edofMat,ones(1,24))',24*24*nele,1);                   % 列索引：


% PREPARE FILTER
iH = ones(nele*(2*(ceil(rmin)-1)+1)^2,1);                                           % 记录当前单元i编号
jH = ones(size(iH));                                                                % 记录邻域单元j编号                                                        
sH = zeros(size(iH));                                                               % 记录单元对权重
k = 0;                                                                              % 计数器：记录单元对数量
for k1 = 1:nelz                                                                     % 遍历三维网格中所有单元
    for i1 = 1:nelx
        for j1 = 1:nely
            e1 = (k1-1)*nelx*nely + (i1-1)*nely+j1;                                 % （k1,i1,j1）位置单元编号：三维坐标转为一维索引编号
            for k2 = max(k1-(ceil(rmin)-1),1):min(k1+(ceil(rmin)-1),nelz)           % 遍历邻域单元：（k2,i2,j2）
                for i2 = max(i1-(ceil(rmin)-1),1):min(i1+(ceil(rmin)-1),nelx)
                    for j2 = max(j1-(ceil(rmin)-1),1):min(j1+(ceil(rmin)-1),nely)
                        e2 = (k2-1)*nelx*nely + (i2-1)*nely+j2;                     % 邻域单元e2的编号
                        k = k+1;                                                    % 已处理的单元对（i,j)数量
                        iH(k) = e1;
                        jH(k) = e2;
                        sH(k) = max(0,rmin-sqrt((i1-i2)^2+(j1-j2)^2+(k1-k2)^2));    % 权重
                    end
                end
            end
        end
    end
end
H = sparse(iH,jH,sH);                                                   % 权重矩阵
Hs = sum(H,2);                                                          % 对H按列求和


% 迭代初始化
x = repmat(volfrac,[nely,nelx,nelz]);           % 初始设计变量矩阵：
xPhys = x;                                      % 物理密度
loop = 0;                                       % 迭代循环次数
change = 1;                                     % 密度变化量


% 开始迭代
while change > tolx && loop < maxloop
    loop = loop+1;
    % 有限元分析
    sK = reshape(KE(:)*(Emin+xPhys(:)'.^penal*(E0-Emin)),24*24*nele,1);
    K = sparse(iK,jK,sK); K = (K+K')/2;
    U(freedofs,:) = K(freedofs,freedofs)\F(freedofs,:);                      % 位移矩阵
    
    % 目标函数与灵敏度分析
    ce = reshape(sum((U(edofMat)*KE).*U(edofMat),2),[nely,nelx,nelz]);       % 最小柔度
    c = sum(sum(sum((Emin+xPhys.^penal*(E0-Emin)).*ce)));
    dc = -penal*(E0-Emin)*xPhys.^(penal-1).*ce;
    dv = ones(nely,nelx,nelz);

    % FILTERING AND MODIFICATION OF SENSITIVITIES
    dc(:) = H*(dc(:)./Hs);  
    dv(:) = H*(dv(:)./Hs);
    % OPTIMALITY CRITERIA UPDATE
    l1 = 0; l2 = 1e9; move = 0.2;
    while (l2-l1)/(l1+l2) > 1e-3
        lmid = 0.5*(l2+l1);
        xnew = max(0,max(x-move,min(1,min(x+move,x.*sqrt(-dc./dv/lmid)))));
        xPhys(:) = (H*xnew(:))./Hs;
        if sum(xPhys(:)) > volfrac*nele, l1 = lmid; else l2 = lmid; end
    end
    change = max(abs(xnew(:)-x(:)));
    x = xnew;

    % 打印结果
    fprintf(' It.:%5i Obj.:%11.4f Vol.:%7.3f ch.:%7.3f\n',loop,c,mean(xPhys(:)),change);

    % 可视化密度结果
    if displayflag, clf; display_3D(xPhys); end %#ok<UNRCH>
end
clf; display_3D(xPhys);
end


% === 生成单元刚度矩阵 ===
function [KE] = lk_H8(nu)

A = [32 6 -8 6 -6 4 3 -6 -10 3 -3 -3 -4 -8;
    -48 0 0 -24 24 0 0 0 12 -12 0 12 12 12];
k = 1/144*A'*[1; nu];

K1 = [k(1) k(2) k(2) k(3) k(5) k(5);            % 子矩阵：不同部分的刚度矩阵
    k(2) k(1) k(2) k(4) k(6) k(7);
    k(2) k(2) k(1) k(4) k(7) k(6);
    k(3) k(4) k(4) k(1) k(8) k(8);
    k(5) k(6) k(7) k(8) k(1) k(2);
    k(5) k(7) k(6) k(8) k(2) k(1)];
K2 = [k(9)  k(8)  k(12) k(6)  k(4)  k(7);
    k(8)  k(9)  k(12) k(5)  k(3)  k(5);
    k(10) k(10) k(13) k(7)  k(4)  k(6);
    k(6)  k(5)  k(11) k(9)  k(2)  k(10);
    k(4)  k(3)  k(5)  k(2)  k(9)  k(12)
    k(11) k(4)  k(6)  k(12) k(10) k(13)];
K3 = [k(6)  k(7)  k(4)  k(9)  k(12) k(8);
    k(7)  k(6)  k(4)  k(10) k(13) k(10);
    k(5)  k(5)  k(3)  k(8)  k(12) k(9);
    k(9)  k(10) k(2)  k(6)  k(11) k(5);
    k(12) k(13) k(10) k(11) k(6)  k(4);
    k(2)  k(12) k(9)  k(4)  k(5)  k(3)];
K4 = [k(14) k(11) k(11) k(13) k(10) k(10);
    k(11) k(14) k(11) k(12) k(9)  k(8);
    k(11) k(11) k(14) k(12) k(8)  k(9);
    k(13) k(12) k(12) k(14) k(7)  k(7);
    k(10) k(9)  k(8)  k(7)  k(14) k(11);
    k(10) k(8)  k(9)  k(7)  k(11) k(14)];
K5 = [k(1) k(2)  k(8)  k(3) k(5)  k(4);
    k(2) k(1)  k(8)  k(4) k(6)  k(11);
    k(8) k(8)  k(1)  k(5) k(11) k(6);
    k(3) k(4)  k(5)  k(1) k(8)  k(2);
    k(5) k(6)  k(11) k(8) k(1)  k(8);
    k(4) k(11) k(6)  k(2) k(8)  k(1)];
K6 = [k(14) k(11) k(7)  k(13) k(10) k(12);
    k(11) k(14) k(7)  k(12) k(9)  k(2);
    k(7)  k(7)  k(14) k(10) k(2)  k(9);
    k(13) k(12) k(10) k(14) k(7)  k(11);
    k(10) k(9)  k(2)  k(7)  k(14) k(7);
    k(12) k(2)  k(9)  k(11) k(7)  k(14)];
KE = 1/((nu+1)*(1-2*nu))*...                % 单元刚度矩阵：子矩阵按照一定方式拼接得到，24*24（矩阵中每个元素表示该单元在特定自由度之间的刚度关系）
    [ K1  K2  K3  K4;
    K2'  K5  K6  K3';
    K3' K6  K5' K2';
    K4  K3  K2  K1'];
end

% === 可视化三维拓扑优化结果 ===
function display_3D(rho)
[nely,nelx,nelz] = size(rho);      % 获取密度矩阵大小
hx = 1; hy = 1; hz = 1;            % 用户自定义单元尺寸
face = [1 2 3 4; 2 6 7 3; 4 3 7 8; 1 5 8 4; 1 2 6 5; 5 6 7 8];      % 单元六个面：每个面有4个顶点
set(gcf,'Name','ISO display','NumberTitle','off');                  % 设置图形窗口属性
for k = 1:nelz                                                      % 遍历所有单元
    z = (k-1)*hz;                                                   % z轴起始位置：从0开始
    for i = 1:nelx
        x = (i-1)*hx;                                               % x轴起始位置
        for j = 1:nely
            y = nely*hy - (j-1)*hy;                                 % y轴起始位置：从顶部开始
            if (rho(j,i,k) > 0.5)  % 用户自定义密度阈值
                vert = [x y z; x y-hx z; x+hx y-hx z; x+hx y z; x y z+hx;x y-hx z+hx; x+hx y-hx z+hx;x+hx y z+hx];                          % 生成当前单元8个顶点坐标
                vert(:,[2 3]) = vert(:,[3 2]); vert(:,2,:) = -vert(:,2,:);                                                                  % 调整顶点坐标的坐标轴顺序和方向：x向右，y向前，z向上
                patch('Faces',face,'Vertices',vert,'FaceColor',[0.2+0.8*(1-rho(j,i,k)),0.2+0.8*(1-rho(j,i,k)),0.2+0.8*(1-rho(j,i,k))]);     % 绘制六面体：根据密度值动态设置灰度深浅
                hold on;
            end
        end
    end
end
axis equal; axis tight; axis off; box on; view([30,30]); pause(1e-6);       % 设置三维视图参数
end
% =========================================================================
% === This code was written by K Liu and A Tovar, Dept. of Mechanical   ===
% === Engineering, Indiana University-Purdue University Indianapolis,   ===
% === Indiana, United States of America                                 ===
% === ----------------------------------------------------------------- ===
% === Please send your suggestions and comments to: kailiu@iupui.edu    ===
% === ----------------------------------------------------------------- ===
% === The code is intended for educational purposes, and the details    ===
% === and extensions can be found in the paper:                         ===
% === K. Liu and A. Tovar, "An efficient 3D topology optimization code  ===
% === written in Matlab", Struct Multidisc Optim, 50(6): 1175-1196, 2014, =
% === doi:10.1007/s00158-014-1107-x                                     ===
% === ----------------------------------------------------------------- ===
% === The code as well as an uncorrected version of the paper can be    ===
% === downloaded from the website: http://www.top3dapp.com/             ===
% === ----------------------------------------------------------------- ===
% === Disclaimer:                                                       ===
% === The authors reserves all rights for the program.                  ===
% === The code may be distributed and used for educational purposes.    ===
% === The authors do not guarantee that the code is free from errors, a